function con = bcfuncFreeTf(s0, s1, el0, el1, mu,cexh,initialThrust)
%s0 = [r(0),theta(0).phi(0),vr(0),vtheta(0),vphi(0),Lamr(0),Lamtheta(0),Lamphi(0),LamVr(0),LamVtheta(0),LamVphi(0),LamTa(0)]
%s1 = same with t_f
ta0 = pi/2;
rtp0 = elements2rtp(el0, ta0, mu);
con(1) = s0(1) - rtp0(1);
con(2) = s0(2) - rtp0(2);
con(3) = s0(3)-rtp0(3);
con(4) = s0(4) - rtp0(4);
con(5) = s0(5) - rtp0(5);
con(6) = s0(6) - rtp0(6);
con(7) = s0(7) - initialThrust;
ta1 = atan2(sin(s1(3)-el1.raan), cos(el1.i)*cos(s1(3)-el1.raan)) - el1.aop;
rtp1 = elements2rtp(el1, ta1, mu);
con(8)  = s1(1) - rtp1(1);
con(9) = s1(2) - rtp1(2);
con(10) = s1(10);
con(11) = s1(4) - rtp1(4);
con(12) = s1(5) - rtp1(5);
con(13) = s1(6) - rtp1(6);
con(14) = s1(14) - 1;
alpha = atan2(-s1(11),-s1(13));
beta = atan2(s1(12), -(s1(11)*sin(alpha) + s1(13)*cos(alpha)));
con(15) = s1(8)*s1(4)+...
    s1(9)*(s1(5)/s1(1))+...
    s1(10)*(s1(6)/(s1(1)*sin(s1(2))))+...
    s1(11)*(s1(5)^2/s1(1) + s1(6)^2/s1(1) - mu/s1(1)^2 + s1(7)*sin(alpha)*cos(beta))+...
    s1(12)*(-s1(4)*s1(5)/s1(1) + s1(6)^2/(s1(1)*tan(s1(2))) - s1(7)*sin(beta))+...
    s1(13)*(-s1(4)*s1(6)/s1(1) - s1(5)*s1(6)/(s1(1)*tan(s1(2))) + s1(7)*cos(alpha)*cos(beta))+...
    s1(14)*(s1(7)^2/cexh);
con = con';
end